


% Properties 
colWidCM = [85, 114, 174]/10; % in mm
% Axes Size
axAspectRatio = 10.5/13;
axWidth = 2.2; % From playing around in illustrator
axHeight = 6.2/3; % Playing around in illustrator
% Markup properties (smaller for insets!)
fontSize = 7;
lineWidth = 0.5;
markerSize = 4;

% Inherited from extracellular plots
%figure(1); clf; hold on;
%set(gcf,'units','centimeters','position',[47 30 axWidth axHeight]);
%set(gca,'units','centimeters','position',[0.05 0.05 axWidth*0.92 axHeight*0.95]); % Note difference for equal axes!

insetWidth = axWidth/2; 
insetHeight = axHeight/3;

dpath = '/Users/landauland/Documents/Research/SabatiniLab/data/manuscriptPreparation/compileData';
fig5results = load(fullfile(dpath,'figure5results.mat'));
res = fig5results.res;

combData = res.combData;
allCa = res.allCa;
allV = res.allV;
patch_amp = res.patch_amp;
mzp = res.mzp;

cdCaData = combData.allCaTraces(:,combData.goodDat) .* combData.dFoverFCa(combData.goodDat);
cdVoData = combData.allVTraces(:,combData.goodDat) .* combData.dFoverFV(combData.goodDat);
idxDistance = combData.distance(combData.goodDat) >= 75;

% Panel A Calcium/Voltage Data
pACa = cdCaData(:,end-4:end);
pAVo = cdVoData(:,end-4:end);
zIdxCa = 10; % Zero Idx (timing of AP)
zIdxVo = 21;



% Calcium Data --- 
tca = (0:16)*50;  % timing = 5
cUseIdx = 3:14;
useTca = tca(cUseIdx) - tca(cUseIdx(1));
caTraces = cell2mat(cellfun(@(c) mean(c,2), {allCa(:).spikeMat}, 'uni', 0));
bcaTraces = caTraces(cUseIdx,:) - mean(caTraces(1:3,:),1);

% Voltage Data -- 
tva = 0:65; % timing = 16
vaTraces = cell2mat(cellfun(@(c) mean(c,2), {allV(:).spikeMat}, 'uni', 0));
bvaTraces = vaTraces - mean(vaTraces(1:10,:),1);
bvaValue = mean(vaTraces(1:10,:),1);

caXScale = tva(end)/useTca(end); % to make them the same width
caYScale = 3.5;
nFOV = size(caTraces,2);

% Generate Figure Insets
axWidth = 5.5; % Based on playing around in illustrator
axHeight = 6.25; % Also based on playing around in illustrator 

doPadMzp = true;
if doPadMzp
    widthReHeight = axWidth / axHeight; 
    mzpSize = size(mzp);
    ogMzpScale = mzpSize(2)/mzpSize(1);
    if ogMzpScale > widthReHeight, error('have to code it differently'); end
    padScale = widthReHeight / ogMzpScale;
    padUnits = floor((padScale * mzpSize(2) - mzpSize(2))/2);
    mzpPad = [0*ones(mzpSize(1),2*padUnits,3), mzp, 0*ones(mzpSize(1),0*padUnits,3)];
else
    mzpPad = mzp;
end

% Calcium Trial CleanUp
idxCTrials = {...
    [1 0 1 1 1 1 1 0 0 1],...
    [1 1 1 0 1 1 1 1 1 1 1],...
    [1 1 1 0 0 0 1 1 1 1],...
    [0 1 1 1 1 1 1],...
    [1 1 1 1 1 1 1 1 1 0 1]};

invertMZP = true;
if invertMZP
    mzpPad = 255 - mzpPad;
end

figure(1); clf; hold on;
set(gcf,'units','centimeters','position',[37 17 axWidth(1) axHeight(1)]);
set(gca,'units','centimeters','position',[0 0 axWidth(1) axHeight(1)]); % Note difference for equal axes!
imagesc(mzpPad);
axis equal
set(gca,'ydir','reverse')

% Plot insets over MZP
pointAxes = [... 
    163.5 333;
    157.5 285;
    184.5 294.6;
    128 250;
    41 153.5];
if doPadMzp
    pointAxes = pointAxes + [2*padUnits 0];
end
for i = 1:nFOV
    plot(pointAxes(i,1),pointAxes(i,2),'marker','.','markersize',markerSize*2,'color','r');
end
% Inset of traces!!!
xTrans = @(x,xShift) 1.7*x + xShift;
yTrans = @(v,yShift) -18*v + yShift;
insetShift = [...
    -127 476 47;
    -127 239.9375 52;
    -127 344.4463 40;
    -127 104.1946 30;
    -127 22.42 24]; 

if doPadMzp
    insetShift = insetShift + [2*padUnits 0 0];
end

allVoltageTrials = cell(1,size(insetShift,1),1);
allCalciumTrials = cell(1,size(insetShift,1),1);
for i = 1:size(insetShift,1)
    cVoltageTrials = allV(i).spikeMat;
    bVoltageTrials = cVoltageTrials - mean(cVoltageTrials(1:10,:),1);
    sVoltageTrials = smoothsmooth(bVoltageTrials,2);
    allVoltageTrials{i} = sVoltageTrials;

    cCalciumTrials = allCa(i).spikeMat(:,logical(idxCTrials{i}));
    bCalciumTrials = cCalciumTrials - mean(cCalciumTrials(cUseIdx(1:3),:),1);
    sCalciumTrials = smoothsmooth(bCalciumTrials,1);
    allCalciumTrials{i} = sCalciumTrials;
end
% plot shaded background to show traces
xExtension = 5;
yExtension = [12 8];
spaceInBetween = zeros(5,2);
for i = 1:size(insetShift,1)
    xLimitBox = insetShift(i,1) + [-xExtension, xTrans(tva(end),0)+xExtension];
    yLimitBox = [(insetShift(i,2)-yExtension(1)+yTrans(max(allCalciumTrials{i},[],'all')/caYScale,0))*[1,1] (insetShift(i,2) + insetShift(i,3) + yTrans(min(allVoltageTrials{i}(1:40,:),[],'all'),0) + yExtension(2))*[1,1]];
    patch([xLimitBox fliplr(xLimitBox)],yLimitBox,0.9*[1,1,1],'EdgeColor','none');
    spaceInBetween(i,:) = yLimitBox([1,3]);
    % and plot connections from ROI to traces
    line([pointAxes(i,1) insetShift(i,1)+xTrans(tva(end),0)+xExtension],[pointAxes(i,2) mean(yLimitBox)],'color','r','linewidth',lineWidth/2);
    line([pointAxes(i,1) insetShift(i,1)+xTrans(tva(end),0)+xExtension],[pointAxes(i,2) mean(yLimitBox)],'color','b','linewidth',lineWidth/2,'linestyle',':');
end
% Used to space insets properly
sib = sort(spaceInBetween);
sum(sib(2:5,1)-sib(1:4,2));
caColTrials = [0.7 0.7 1];
caCol = 'b';
plotTrials = true;
for i = 1:size(insetShift,1)%nFOV
    if plotTrials
        plot(xTrans(tva,insetShift(i,1)),yTrans(allVoltageTrials{i},insetShift(i,2)+insetShift(i,3)),'color',[1 0.7 0.7],'linewidth',lineWidth/2)
        plot(xTrans(tva,insetShift(i,1)),mean(yTrans(allVoltageTrials{i},insetShift(i,2)+insetShift(i,3)),2),'color','r','linewidth',lineWidth)
        plot(xTrans(useTca*caXScale,insetShift(i,1)),yTrans(allCalciumTrials{i}(cUseIdx,:)/caYScale,insetShift(i,2)),'color',caColTrials,'linewidth',lineWidth/3)
        plot(xTrans(useTca*caXScale,insetShift(i,1)),yTrans(mean(allCalciumTrials{i}(cUseIdx,:)/caYScale,2),insetShift(i,2)),'color',caCol,'linewidth',lineWidth)
    else
        plot(xTrans(tva,insetShift(i,1)),yTrans(bvaTraces(:,i),insetShift(i,2)),'color','r','linewidth',lineWidth/2)
        plot(xTrans(useTca*caXScale,insetShift(i,1)),yTrans(bcaTraces(:,i)/caYScale,insetShift(i,2)+insetShift(i,3)),'color',caCol,'linewidth',lineWidth/2)
    end
end
% Create Legend & Scale Bar
vScalePos = [245 175];
if doPadMzp
    vScalePos = vScalePos + [2*padUnits 0];
end
line(xTrans([0 20],vScalePos(1)),vScalePos([2,2]),'color','r','lineWidth',lineWidth);
line(vScalePos([1,1]),yTrans([0,1],vScalePos(2)),'color','r','lineWidth',lineWidth);
text(vScalePos(1),vScalePos(2),'20 ms','color','r','fontSize',fontSize,'HorizontalAlignment','left','VerticalAlignment','top');
text(vScalePos(1),vScalePos(2),'10% \DeltaF/F','color','r','fontSize',fontSize,'HorizontalAlignment','left','VerticalAlignment','bottom','Rotation',90);
text(vScalePos(1)-22,vScalePos(2),'Voltage','FontSize',fontSize,'Color','r','HorizontalAlignment','Left','VerticalAlignment','Bottom','Rotation',90);
% Moving scale bar
cScalePos = [140 175];
if doPadMzp
   cScalePos = cScalePos + [2*padUnits 0];
end
line(xTrans([0 200]*caXScale,cScalePos(1)),cScalePos([2,2]),'color',caCol,'lineWidth',lineWidth);
baseline2MaxBiggest = [0 8.8]; % This is the baseline to max value for ROI 2 (which is the biggest one)
actualDffBiggest = pACa(zIdxCa,1); % The 
line(cScalePos([1,1]),yTrans([0,8.8]/caYScale*0.3/actualDffBiggest,cScalePos(2)),'color',caCol,'lineWidth',lineWidth); %
text(cScalePos(1),cScalePos(2),'200 ms','color',caCol,'fontSize',fontSize,'HorizontalAlignment','left','VerticalAlignment','top');
text(cScalePos(1),cScalePos(2),'30% \DeltaF/F','color',caCol,'fontSize',fontSize,'HorizontalAlignment','left','VerticalAlignment','bottom','Rotation',90);
text(cScalePos(1)-22,cScalePos(2),'Calcium','FontSize',fontSize,'Color',caCol,'HorizontalAlignment','Left','VerticalAlignment','Bottom','Rotation',90);
tvmScale = 10/1000; % 10 ms / 1000 samples
vmScale = 20/0.4; % 20 mV / 0.4 value
tvm = (0:6500) * tvmScale;
useIdx = 500:6000;
useTvm = tvm(useIdx);
wcXScale = 1.5;
wcYScale = 10;
whiteCol = 0*[1,1,1];
%Repeat the same nBack and nFront of the camera (in ms)
spikeT_patch = spikefindhyst(patch_amp, 0.5, -0.6);
daq_samplerate = 1e5;
nBack = 0.020*daq_samplerate; nFront = 0.045*daq_samplerate;
nSpike_patch = length(spikeT_patch);
spikeMat_patch = zeros(nBack + nFront + 1, nSpike_patch);
for k = 1:nSpike_patch
    spikeMat_patch(:,k) = patch_amp((spikeT_patch(k)-nBack):(spikeT_patch(k)+nFront));
    spikeMat_patch(:,k) = spikeMat_patch(:,k) - prctile(spikeMat_patch(:,k), 10);
end
avgSpikeWaveform = smoothsmooth(mean(spikeMat_patch,2),20)*vmScale;
plot(xTrans(useTvm,205+2*padUnits*doPadMzp),yTrans(avgSpikeWaveform(useIdx)/wcYScale,500),'color',whiteCol,'linewidth',lineWidth/2);
wScalePos = [268 450];
if doPadMzp
    wScalePos = wScalePos + [2*padUnits 0];
end
line(xTrans([0 20],wScalePos(1)),wScalePos([2,2]),'color',whiteCol,'lineWidth',lineWidth);
line(wScalePos([1,1]),yTrans([0,20]/wcYScale,wScalePos(2)),'color',whiteCol,'lineWidth',lineWidth);
text(wScalePos(1),wScalePos(2),'20 ms','color',whiteCol,'fontSize',fontSize,'HorizontalAlignment','left','VerticalAlignment','top');
text(wScalePos(1),wScalePos(2),'20 mV','color',whiteCol,'fontSize',fontSize,'HorizontalAlignment','left','VerticalAlignment','bottom','Rotation',90);
set(gca,'visible','off');



dtV = 1;
dtCa = 50;  % set to 1 to have equal axes.  Actual time interval is 50 ms
vTau = (-20:45)*dtV;
CaTau = (-9:15)*dtCa;
cdCaDistData = cdCaData(:,idxDistance);
cdVoDistData = cdVoData(:,idxDistance);
caPeak = cdCaDistData(zIdxCa,:); caPeak(caPeak<0.001)=0.001;
voPeak = cdVoDistData(zIdxVo,:);
% Convert range of selection to the range of DF/F Cohen lab measured
% (this is just for visualization)
caScaleRange = 1.75;
idxSilent = caPeak <= 0.04*caScaleRange;
idxActive = caPeak >= 0.1*caScaleRange  &  caPeak <= 0.3*caScaleRange;
xLim = [min(CaTau) max(CaTau)];
yLim = [-0.05 0.2825];
axWidth = 2.2; % From playing around in illustrator
axHeight = 6.2/3; % Playing around in illustrator


% Figure 5B
figure(2); clf; hold on;
set(gcf,'units','centimeters','position',[47 30 axWidth axHeight]);
set(gca,'units','centimeters','position',[0.1 0.05 axWidth*0.965 axHeight*0.965]); % Note difference for equal axes!
shadedErrorBar(CaTau, mean(cdCaDistData(:,idxActive),2),std(cdCaDistData(:,idxActive),1,2)/sqrt(sum(idxActive)),{'color','k','linewidth',lineWidth},0);
shadedErrorBar(CaTau, mean(cdCaDistData(:,idxSilent),2),std(cdCaDistData(:,idxActive),1,2)/sqrt(sum(idxSilent)),{'color','b','linewidth',lineWidth},0);
xlim(xLim);
ylim(yLim);
% Mark AP Time
line(-[100 100],[-0.045 -0.025],'color','k','linewidth',lineWidth);
% Scale Bar
scalePosition = [320 0.09 150 0.05];
line(scalePosition(1)+[0 scalePosition(3)],scalePosition(2)*[1,1],'linewidth',lineWidth,'color','k');
line(scalePosition(1)*[1,1],scalePosition(2)+[0 scalePosition(4)],'linewidth',lineWidth,'color','k');
text(scalePosition(1),scalePosition(2),sprintf('%d ms',scalePosition(3)),'Fontsize',fontSize,'HorizontalAlignment','Left','VerticalAlignment','Top');
text(scalePosition(1),scalePosition(2),[sprintf('%d%%',100*scalePosition(4)),' \DeltaF/F'],'Fontsize',fontSize,'HorizontalAlignment','Left','VerticalAlignment','Bottom','Rotation',90);
set(gca,'visible','off');


% Figure 5C
xLim = [min(vTau) max(vTau)];
yLim = [-0.05*0.4/0.2825 0.4];
figure(3); clf; hold on;
set(gcf,'units','centimeters','position',[47 23.5 axWidth axHeight]);
set(gca,'units','centimeters','position',[0.1 0.05 axWidth*0.965 axHeight*0.965]); % Note difference for equal axes!
shadedErrorBar(vTau, mean(cdVoDistData(:,idxActive),2),std(cdVoDistData(:,idxActive),1,2)/sqrt(sum(idxActive)),{'color','k','linewidth',lineWidth},0);
shadedErrorBar(vTau, mean(cdVoDistData(:,idxSilent),2),std(cdVoDistData(:,idxActive),1,2)/sqrt(sum(idxSilent)),{'color','r','linewidth',lineWidth},0);
xlim(xLim);
ylim(yLim);
% Mark AP Time
line(-[7 7],[-0.045 -0.025]*0.4/0.2825,'color','k','linewidth',lineWidth);
% Scale Bar
scalePosition = [20 0.135 5 0.05];
line(scalePosition(1)+[0 scalePosition(3)],scalePosition(2)*[1,1],'linewidth',lineWidth,'color','k');
line(scalePosition(1)*[1,1],scalePosition(2)+[0 scalePosition(4)],'linewidth',lineWidth,'color','k');
text(scalePosition(1),scalePosition(2),sprintf('%d ms',scalePosition(3)),'Fontsize',fontSize,'HorizontalAlignment','Left','VerticalAlignment','Top');
text(scalePosition(1),scalePosition(2),[sprintf('%d%%',100*scalePosition(4)),' \DeltaF/F'],'Fontsize',fontSize,'HorizontalAlignment','Left','VerticalAlignment','Bottom','Rotation',90);
set(gca,'visible','off');



% Figure 5D
setBackX = [-0.05 -0.1 -0.07 -0.42]/2*0.9;
setBackY = [-0.015 -0.028 -0.032 -0.155]/0.45*0.7;
yLim = [0 0.95];
figure(14); clf; hold on;
set(gcf,'units','centimeters','position',[47 17 axWidth axHeight]);
set(gca,'units','centimeters','position',[0.65 0.6 axWidth*0.71 axHeight*0.7]); % Note difference for equal axes!
% Red/Blue Patch
patchHeight = diff(yLim)/9;
for i = yLim(1):patchHeight:yLim(2)-patchHeight
    ccol = 'b'; if rem(int8((i-yLim(1))/patchHeight),2)==0, ccol = 'r'; end
    patch([0 0.04 0.04 0]*caScaleRange,i+[0 0 patchHeight patchHeight],ccol,'FaceColor',ccol,'FaceAlpha',0.2,'EdgeColor','none');
end
patch([0.1 0.3 0.3 0.1]*caScaleRange,[0 0 0.95 0.95],'k','FaceColor','k','FaceAlpha',0.2,'EdgeColor','none');
plot(caPeak, voPeak,'marker','.','color','k','linestyle','none','markersize',markerSize);
xlim([setBackY(2) 0.75]);
ylim([setBackX(2) 0.95]);
line([0 0.75],setBackX(1)*[1,1],'color','k','linewidth',lineWidth);
line(setBackY(1)*[1,1],[0 0.95],'color','k','linewidth',lineWidth);
xticks = 0:0.3:0.6;
for xt = xticks
    line(xt*[1,1],setBackX([1 2]),'color','k','linewidth',lineWidth);
    text(xt,setBackX(3),num2str(xt),'Fontsize',fontSize,'HorizontalAlignment','Center','VerticalAlignment','Top');
end
text(0.7/2, setBackX(4), '\DeltaCa_{AP} (\DeltaF/F)','Fontsize',fontSize,'HorizontalAlignment','Center','VerticalAlignment','Top');
yticks = 0:0.4:0.8;
for yt = yticks
    line(setBackY([1 2]),yt*[1,1],'color','k','linewidth',lineWidth);
    text(setBackY(3),yt,num2str(yt),'Fontsize',fontSize,'HorizontalAlignment','Right','VerticalAlignment','Middle');
end
text(setBackY(4), 0.95/2, '\DeltaV (\DeltaF/F)','Fontsize',fontSize,'HorizontalAlignment','Center','VerticalAlignment','Bottom','Rotation',90);
set(gca,'visible','off');








%% -- quantify time course of voltage traces (for resubmission) -- 

[pkValue,idxPeak] = max(cdVoDistData,[],1);
mdFiltVo = medfilt1(cdVoDistData,3,[],1);
halfWidthSamples = zeros(length(pkValue),1);
for roi = 1:length(pkValue)
    leftSamples = find(mdFiltVo(idxPeak(roi):-1:1,roi)<=pkValue(roi)/2,1);
    rightSamples = find(mdFiltVo(idxPeak(roi):end,roi)<=pkValue(roi)/2,1);
    halfWidthSamples(roi) = leftSamples + rightSamples - 2; 
end
halfWidthTime = halfWidthSamples * mean(diff(vTau)); % (dtau is 1)

fprintf('Active HW: %.1f +/- %.1f ms\n', mean(halfWidthTime(idxActive)), std(halfWidthTime(idxActive)));
fprintf('Silent HW: %.1f +/- %.1f ms\n', mean(halfWidthTime(idxSilent)), std(halfWidthTime(idxSilent)));


